cx = [-1 1 1 -1];
cy = [1 1 -1 -1];

bx = [-1 1 1 -1];
by = [1 2 0 -1];

[x, y] = generateBoundaryNode(cx, cy, 100);
[dx, dy] = generateBoundaryNode(bx, by, 100);

domain2D = createDomain(x, y);

dim = size(domain2D.x, 2);

b = zeros(2*dim, 1);
b(1:2:2*dim) = dx' - x';
b(2:2:2*dim) = dy' - y';

% rx = rand(100, 1) * 2 - 1;
% ry = rand(100, 1) * 2 - 1;

rx = [-0.5; -0.5; 0.5; 0.5];
ry = [-0.5; 0.5; 0.5; -0.5];

A = assemblyMatrix(domain2D);

% b = assemblyBoundaryVector(domain2D);

theta = A \ b;

u = evaluateInteriorPoint(rx, ry, domain2D, theta);

plotDeformation(domain2D, rx, ry, b(1:2:2*dim), b(2:2:2*dim), u(:,1), u(:,2));

% dx = rand(1, 1);
% dy = rand(1, 1);
% 
% rc = rand(1000*9, 1) * 100;
% strain = reshape(rc, 3, 3, 1000);
% 
% plotDeformation(domain2D, rx, ry, zeros(dim, 1), zeros(dim, 1), dx, dy);